library(parallel)
library(doMC)
library(DHARMa)
library(broom.mixed)
library(MuMIn)
library(mgcv)
library(bbmle)
library(lme4)
library(merTools)
library(magrittr)
library(gridExtra)
library(ggplot2)
library(tidyr)
library(purrr)
library(plyr)
library(dplyr)1 - Variables
table 1: Variáveis utilizadas na descrição estatística
| code | description | class | range |
|---|---|---|---|
| n_nRef | number of SADs not refuted by the KS bootstrap test with critical p-value 0.05 | continuous | (0 ; 100) |
| diffS_mean | diffS = (S_MN - S_obs) / S_obs, diffS_mean = mean(diffS); S_MN = species richness estimated by MN | continuous | (-0.977 ; 4.78) |
| U_med | average of 10 replicates of the speciation rate estimated by MNEE | continuous | (8.860e-5 ; 4.221e-2) |
| p | proportion of tree cover | continuous | (0.013 ; 0.961) |
| MN | Neutral Model (EE = spatial explicit; EI = spatial implict) | category | 2 levels |
| d | mean dispersal distance (meters) | continuous | (0.456 ; 19.183) |
| k | proportion of propagules until the nearest neighborhood seq(0.99 : 0.05) | category | 20 levels |
| d_Lplot | mean dispersal distance / side of the sample area | continuous | (0.02 ; 0.192) |
| S_obs | observed species richness | integer | (26 ; 458) |
| Ntotal | number of individuals in the sample area | integer | (649 ; 12105) |
| SiteCode | sample area code | category | 103 levels |
Figure 1 Possible Predictor Variables
2 - Number of unrefuted SADs
Figura 2.1 número de SADs não refutadas ~ p * k * MN. A linha azul é uma estimativa baseada em ‘loess’.
Figura 2.2 número de SADs não refutadas ~ (d * MN|Site ~ p). Site está ordenado pelo valor de p.
2.1 GLMM binomial
Modelo cheio 1:
linear term
- ~ (p + p^2) * (d + d^2) * MN; if : var_dispersao = contiguous
- ~ ( p + p^2) * k * MN; if : var_dispersao = category
random term
- 1|SiteCode
- MN|SiteCode
- (var_dispersao + var_dispersao^2)*MN|SiteCode, if : var = dispersao : contiguous
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",3)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio) <- NULL
l_nRef__modeloCheio <- do.call(c,l_nRef__modeloCheio)
#
f_warning <- function(glmm_object){
# glmm_object <- l_nRef__modeloCheio[[4]]
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.1 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00343548 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | Model failed to converge with max|grad| = 0.00836802 (tol = 0.002, component 1) |
| d_Lplot.z (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 40.1827 (tol = 0.002, component 1) |
| d.z MN|Site | Model failed to converge with max|grad| = 0.00746139 (tol = 0.002, component 1) |
| d.z (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
i <- 1
names(l_nRef__modeloCheio[v_glmerUpdate])[i]
#
## update1:
md_allFit <- allFit(l_nRef__modeloCheio[v_glmerUpdate][[i]],maxfun = 1e5, parallel = 'multicore', ncpus = 3)1) d_Lplot 1|Site
7 optimizer(s) failed
2) d_Lplot MN|Site
7 optimizer(s) failed
3) d 1|Site
7 optimizer(s) failed
4) d MN|Site
7 optimizer(s) failed
Modelo cheio 2:
# dados
df_md <- df_resultados %>%
select(n_nRef,p,SiteCode,MN,p.z,k,d.z,d_Lplot.z)
df_md[,6:8] <- apply(df_md[,6:8],2,as.character)
df_md %<>%
gather(key = var_dispersao, value = value_dispersao, k:d_Lplot.z)
# funcao ajuste dos modelos
f_md <- function(dados){
var_dispersao <- unique(dados$var_dispersao)
if(var_dispersao == "k"){
l_md <- vector("list",2)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"))
dados$value_dispersao <- factor(dados$value_dispersao,levels = unique(dados$value_dispersao)[20:1])
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}else{
dados$value_dispersao <- as.numeric(dados$value_dispersao)
l_md <- vector("list",4)
names(l_md) <- c(paste0(var_dispersao," 1|Site"),
paste0(var_dispersao," MN|Site"),
paste0(var_dispersao," ",var_dispersao,"*MN|Site"),
paste0(var_dispersao,"^2 (",var_dispersao,"+",var_dispersao,"^2)*MN|Site"))
l_md[[1]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(1|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[2]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[3]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * value_dispersao * MN +
(value_dispersao * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
l_md[[4]] <- glmer(cbind(n_nRef,100-n_nRef) ~
(p.z + I(p.z^2)) * (value_dispersao + I(value_dispersao^2)) * MN +
( (value_dispersao + I(value_dispersao^2) ) * MN|SiteCode),
family = "binomial",data=dados,
control=glmerControl(optimizer="bobyqa",optCtrl=list(maxfun=2e9)))
}
return(l_md)
}
registerDoMC(3)
l_nRef__modeloCheio2 <- dlply(df_md,.(var_dispersao),f_md,.parallel = TRUE)
names(l_nRef__modeloCheio2) <- NULL
l_nRef__modeloCheio2 <- do.call(c,l_nRef__modeloCheio2)
#
f_warning <- function(glmm_object){
v_message <- glmm_object@optinfo$conv$lme4$messages %>%
as.character()
if(length(v_message)==0){v_message <- "OK"}
if(length(v_message)>1){v_message <- v_message[1] }
return(v_message)
}
df_auditoria.md <- ldply(l_nRef__modeloCheio2,f_warning,.id="glmer") %>%
rename(warning_message = V1)Table 2.2 Modelos Cheios estimados e avisos de convergência
| glmer | warning_message |
|---|---|
| d_Lplot.z 1|Site | Model failed to converge with max|grad| = 0.00284637 (tol = 0.002, component 1) |
| d_Lplot.z MN|Site | OK |
| d_Lplot.z d_Lplot.z*MN|Site | OK |
| d_Lplot.z^2 (d_Lplot.z+d_Lplot.z^2)*MN|Site | OK |
| d.z 1|Site | Model failed to converge with max|grad| = 0.00369235 (tol = 0.002, component 1) |
| d.z MN|Site | OK |
| d.z d.z*MN|Site | OK |
| d.z^2 (d.z+d.z^2)*MN|Site | OK |
| k 1|Site | OK |
| k MN|Site | OK |
allFit
d_Lplot.z 1|Site
7 optimizer(s) failed
d.z 1|Site
7 optimizer(s) failed
Comparação de Modelos Cheios:
Tabela 2.3 Comparação baseada em AICc dos modelos cheios
| GLMM | dAICc | df | weight |
|---|---|---|---|
| d^2 (d+d^2)*MN|Site | 0.0 | 39 | 1 |
| d_Lplot^2 (d_Lplot+d_Lplot^2)*MN|Site | 274.6 | 39 | <0.001 |
| d d*MN|Site | 38427.3 | 22 | <0.001 |
| d_Lplot d_Lplot*MN|Site | 38586.5 | 22 | <0.001 |
| k MN|Site | 74141.2 | 123 | <0.001 |
| d MN|Site | 95710.0 | 15 | <0.001 |
| d_Lplot MN|Site | 96279.6 | 15 | <0.001 |
| k 1|Site | 106158.0 | 121 | <0.001 |
Diagnostico do modelo cheio mais plausível 1
Figura 1.2 Resíduos Quantílicos do modelo cheio plausível
Diagnostico do modelo cheio mais plausível 2
